

Numerical Simulation of Turbulent MHD Flows 
Using an Iterative PNS Algorithm 

Hiromasa Kato* and John C. Tannehill 1 
Iowa State University, Ames, IA 50011 
and 

Unmeel B. Mehta 11 

NASA Ames Research Center, Moffett Field, CA 94035 


Abstract 

A new parabolized Navier-Stokes (PNS) algorithm 
has been developed to efficiently compute magne- 
tohydrodynamic (MHD) flows in the low magnetic 
Reynolds number regime. In this regime, the electri- 
cal conductivity is low and the induced magnetic field 
is negligible compared to the applied magnetic field. 
The MHD effects are modeled by introducing source 
terms into the PNS equation which can then be solved 
in a very efficient manner. To account for upstream 
(elliptic) effects, the flowfields are computed using 
multiple streamwise sweeps with an iterated PNS al- 
gorithm. Turbulence has been included by modify- 
ing the Baldwin-Lomax turbulence model to account 
for MHD effects. The new algorithm has been used 
to compute both laminar and turbulent, supersonic, 
MHD flows over flat plates and supersonic viscous 
flows in a rectangular MHD accelerator. The present 
results are in excellent agreement with previous com- 
plete Navier-Stokes calculations. 

Introduction 

Flowfields involving MHD effects have typically been 
computed [1-10] by solving the complete Navier- 
Stokes (N-S) equations for fluid flow in conjunction 
with Maxwell’s equations of electromagnetodynam- 
ics. When chemistry and turbulence effects are also 
included, the computational effort required to solve 
the resulting coupled system of partial differential 
equations is extremely formidable. One possible rem- 
edy to this problem is to use the parabolized Navier- 

* Graduate Research Assistant, Student Member AIAA 

t Manager, Computational Fluid Dynamics Center, and 
Professor, Dept, of AEEM. Fellow AIAA 

* Division Scientist, Associate Fellow, AIAA 

Copyright ©2003 by the American Institute of Aeronau- 
tics and Astronautics, Inc., all rights reserved. 


Stokes (PNS) equations in place of the N-S equa- 
tions. The PNS equations can be used to compute 
three-dimensional, supersonic viscous flowfields in a 
very efficient manner [11]. This efficiency is achieved 
because the equations can be solved using a space- 
marching technique as opposed to the time-marching 
technique that is normally employed for the complete 
N-S equations. 

Recently, the present authors have developed a 
PNS code to solve supersonic MHD flowfields in 
the high magnetic Reynolds number regime [12]. 
This code is based on NASA’s upwind PNS (UPS) 
code which was originally developed by Lawrence et 
al. [13]. The UPS code solves the PNS equations us- 
ing a fully conservative, finite-volume approach in a 
general nonorthogonal coordinate system. The UPS 
code has been extended to permit the computation of 
flowfields with strong upstream influences. In regions 
where strong upstream influences are present, the 
governing equations are solved using multiple sweeps. 
As a result of this approach, a complete flowfield can 
be computed more efficiently (in terms of computer 
time and storage) than with a standard N-S solver 
which marches the entire solution in time. Three it- 
erative PNS algorithms (IPNS, TIPNS, and FBIPNS) 
have been developed. The iterated PNS (IPNS) al- 
gorithm [14] can be applied to flows with moder- 
ate upstream influences and small streamwise sepa- 
rated regions. The time iterated PNS (TIPNS) algo- 
rithm [15] can be used to compute flows with strong 
upstream influences including large streamwise sep- 
arated regions. The forward-backward sweeping it- 
erative PNS (FBIPNS) algorithm [16] was recently 
developed to reduce the number of sweeps required 
for convergence. 

The majority of MHD codes that have been de- 
veloped combine the electromagnetodynamic equa- 
tions with the full Navier-Stokes equations result- 
ing in a complex system of eight scalar equations. 
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= V-(V.?)-V-U + E-J 


These codes can theoretically be used for any mag- 
netic Reynolds number which is defined as Re m = 
cTe/ieKooL where <r e is the electrical conductivity, fi e 
is the magnetic permeability, Voo is the freestream 
velocity, and L is the reference length. However, 
it has been shown that as the magnetic Reynolds 
number is reduced, numerical difficulties are often 
encountered [4]. For many aerospace applications 
the electrical conductivity of the fluid is low and 
hence the magnetic Reynolds number is small. In 
these cases, it makes sense to use the low magnetic 
Reynolds number assumption and reduce the com- 
plexity of the governing equations. In this case, the 
MHD effects can be modeled with the introduction 
of source terms into the fluid flow equations. Severed 
investigators [4,8,17-19] have developed N-S codes 
for the low magnetic Reynolds number regime where 
the induced magnetic field is negligible compared to 
the applied magnetic field. 

In the present study, a new PNS code (based on 
the UPS code) has been developed to compute MHD 
flows in the low magnetic Reynolds number regime. 
The MHD effects are modeled by introducing the ap- 
propriate source terms into the PNS equations. Up- 
stream elliptic effects cam be accounted for by us- 
ing multiple streamwise sweeps with either the IPNS, 
TIPNS, or FBIPNS adgorithms. Turbulence has been 
included by modifying the Baldwin-Lomax turbu- 
lence model [20] to account for MHD effects using 
the approach of Lykoudis [21]. The new code has 
been tested by computing both laminar and turbu- 
lent, supersonic MHD flows over a flat plate. Com- 
parisons have been made with the previous complete 
N-S computations of Dietiker and Hoffmann [18]. In 
addition, the new code has been used to compute the 
supersonic viscous flow inside a rectamgular channel 
designed for MHD experiments [22]. 


( 3 ) 


Ohm's law 

J = (r c (E + VxB) (4) 

where V is the velocity vector, B is the magnetic 
field vector, E is the electric field vector, and J is the 
conduction current density. The flow is assumed to 
be either in chemical equilibrium or in a frozen state. 
The curve fits of Srinivasan et al. [23, 24] are used 
for the thermodynamic and transport properties of 
equilibrium air. 

The governing equations are nondimensionalized 
using the following reference variables. 
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where the superscript * refers to the nondimensional 
quantities. In subsequent sections, the asterisks are 
dropped. 

If the flow variables are assumed to vary in only two 
dimensions (x, y) while the velocity, magnetic, and 
electric fields have components in three dimensions 
(x, y, z), the governing equations can be written in 
the following flux vector form: 


Governing Equations 

The governing equations for a viscous MHD flow with 
a small magnetic Reynolds number are given by [18]: 
Continuity equation 

^ + V - (pV) = 0 (1) 


3U dEi dFi dE v dF„ 
~dt + ~d^ + ~dy = ~dx + ~dy + Smhd 


( 6 ) 


where U is the vector of dependent variables and E, 
and F , are the in viscid flux vectors, and E v and F v 
are the viscous flux vectors. The source term Smhd 
contains all of the MHD effects. The flux vectors are 
given by 


Momentum equation 

d(pV) 


dt 

Energy equation 

djpet) 

dt 


+ V-Lw + p!]=V-f + JxB (2) 
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P e t — x/> ( u2 + v 2 + u;2 ) + c — (12) 

2 . 7 - X 

and 7 can be determined from the curve fits of Srini- 
vasan et al. [23] for an equilibrium air flow or is equal 
to a constant (7) for a frozen or perfect gas flow. The 
nondimensional shear stresses and heat fluxes are de- 
fined in the usual manner [11]. 

The governing equations are transformed into com- 
putational space and written in a generalized coordi- 
nate system (£,77) as 

yU t + E f + F, = ^S (13) 

where 

E = (y)(B*-E.)+^(F i -F.) 

F = (y)(B*-B-)+(^)( F<-F w ) (14) 

and J is the Jacobian of the transformation. 

The governing equations are parabolized by drop- 
ping the time derivative term and the streamwise 
direction (£) viscous flow terms in the flux vectors. 
Equation (13) can then be rewritten as 


+ F„ — 


Smhd 


(15) 


where 

E 
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7 > + (^' f ‘ 


= (y) (E ,• - K) + (&) (F, - F' v ) (16) 


The prime in the preceding equation indicates that 
the streamwise viscous flow terms have been dropped. 

For turbulent flows, the two-layer Baldwin- Lorn ax 
turbulence model [20] has been modified to account 
for MHD effects. Only the expression for turbulent 
viscosity in the inner layer is changed. This modifi- 
cation for MHD flows is due to Lykoudis [9,21]. 


Numerical Method 


The governing PNS equations with MHD source 
terms have been incorporated into NASA’s upwind 
PNS (UPS) code [13]. These equations can be solved 
very efficiently using a single sweep of the flowfield 
for many applications. For cases where upstream 
(elliptic) effects are important, the flowfield can be 
computed using multiple streamwise sweeps with ei- 
ther the IPNS [14], TIPNS [15], or FBIPNS [16] algo- 
rithms. This iterative process is continued until the 
solution is converged. 

For the iterative PNS (IPNS) method, the E vec- 
tor is split using the Vigneron parameter (u;) [25]. 
This parameter does not need to be changed for the 
present low magnetic Reynolds number formulation. 
In the previous high magnetic Reynolds number code 
[12] it was necessary to modify the Vigneron param- 
eter to account for MHD effects. After splitting, the 
E vector can be written as: 
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(18) 


The streamwise derivative of E is then differenced 
using a forward difference for the “elliptic” portion 

(E p ): 


(f ) j+1 = ift( E «- E . r ) + ( E ?«- E W] 

(19) 
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where the subscript (i 4* 1) denotes the spatial index 
(in the f direction) where the solution is currently 
being computed. The vectors E* +1 and Ef +1 are then 
linearized in the following manner: 

E 'hi = E * + (w),. (u *' +1 ~ Ui) 

E? +1 = E ? + (w) i (U<+ 1 - U<) (20) 


The Jacobians can be represented by 


A' 

A p 


dE‘ 

au 

SEP 

au 


( 21 ) 


After substituting the above linearizations into 
Eq. (19), the expression for the streamwise gradient 
of E becomes 



A* 


(AT-Af)( U, +1 -U,) 


+ (E?+2 — Ef) ] (22) 


The final discretized form of the fluid flow equa- 
tions with MHD source terms is obtained by substi- 
tuting Eq. (22) into Eq. (15) along with the linearized 
expression for the flux in the cross flow plane. The 
final expression becomes: 


Test Case 1: Supersonic laminar and 
turbulent flows over a flat plate with 
applied magnetic field 

In this test case, the supersonic, laminar and turbu- 
lent flow over a flat plate with an applied magnetic 
field is computed. This case corresponds to the flat 
plate case computed previously by Dietiker and Hoff- 
mann [18] using the full N-S equations. A strong 
magnetic field is applied normal to the flow as shown 
in Fig. 1. The dimensional flow parameters for this 
test case are: 


Moo 

= 

2.0 

Poo 

= 

1.076 x 10 s N/m : 

Too 

= 

300 K 

Rcoo 

= 

3.75 x 10 6 

7 

= 

1.4 

L 

= 

0.08 m 

<r e 


800 mho/m 


The plate is assumed to be an adiabatic wall and a 
perfect gas flow is assumed. The magnetic Reynolds 
number (based on the length of the plate) is 0.056 and 
cam be considered negligible when compared to one. 
The normal magnetic field component (B y ) ranges 
in value from 0.0 to 1.2 T. The magnitude of the 
magnetic field cam be represented by the parameter 
m which is defined [18] by 



= RHS (23) 

where 

RHS = _J_[ (E f + 2 )*-(E?) fe+1 ] 

and the superscript k-fl denotes the current iteration 
(i.e. sweep) level. 


Numerical Results 

In order to investigate the utility and accuracy of 
the present PNS approach of solving MHD flowfields 
at low magnetic Reynolds numbers, a few basic test 
cases were computed. The supersonic viscous flow in 
these cases is altered by the presence of the magnetic 
and electric fields which are applied to the flow. 


Poo Uoo 


(24) 


and has units of (1/m). For = 1.2 T, m is equal 
to 1.33. 

A highly stretched grid consisting of 50 points in 
the normal direction was used to compute this case. 
The first point off the wall was located at 2 x 10” 7 
m. Initially, the flow was assumed laminar and sev- 
eral values of B y ranging from 0.0 (no magnetic field) 
to 1.2 T were used. The velocity and temperature 
profiles at x = 0.06 m are shown in Figs. 2 and 3 for 
B y = 0.0 T, 1.0 T, and 1.2 T. The velocity profiles 
are compared to the N-S results of Dietiker and Hoff- 
mann in Fig. 2 and show excellent agreement. The 
magnetic field generates a Lorentz force which acts in 
a direction opposite to the flow. Thus, the flow is de- 
celerated as the magnetic field is increased as seen in 
Fig. 2. For B y — 1.2 T the flow is slightly separated. 
The temperature profiles cannot be compared at this 
time since no temperature data is given in Ref. [18]. 

The turbulent flow over the flat plate was then 
computed using the modified Baldwin-Lomax turbu- 
lence model that accounts for MHD effects. The flow 
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was assumed laminar prior to the point (z = 0.04 m) 
where transition from laminar to turbulent flow was 
triggered. Again, several values of B y ranging from 
0.0 to 1.2 T were used in the computations. The tur- 
bulent velocity and temperature profiles at z = 0.06 
m are shown in Figs. 4 and 5 for B y = 0.0, 1.0 T, 
and 1.2 T. The turbulent velocity profiles in Fig. 4 
are in good agreement with the results of Ref. [18]. 
The variation of skin friction coefficient is shown in 
Fig. 6. The present laminar/ turbulent skin friction 
variations are compared with the results of Ref. [18] 
and show' good agreement. The difference in results 
near the transition point may be due to the coarse 
grid and smoothing used in Ref. [18]. 

All of the present laminar computations were per- 
formed using a single sweep of the flowfield except for 
the separated flow case ( B y — 1.2 T). For this case 
as well as for all the turbulent cases, multiple sweeps 
were used to account for upstream effects. 

Test Case 2: Supersonic viscous flow in 
a rectangular MHD accelerator 

In this test case, the supersonic flow in an experi- 
ment ad MHD channel is simulated. This facility is 
currently being built at NASA Ames Research Center 
by D. W. Bogdanoff, C. Park, and U. B. Mehta [22] 
to study criticad technologies related to MHD bypass 
scramjet engines. The channel is about a half meter 
long and contains a nozzle section, a center section, 
and an aiccelerator section. The channel has a uni- 
form width of 2.03 cm. Magnetic and electric fields 
can be imposed upon the flow in the aiccelerator sec- 
tion. A schematic of the MHD accelerator section is 
shown in Fig. 7. 

This test case was previously computed by 
R. W. MacCormack [10] using the full N-S equations 
coupled with the electromagnetodynamic equations. 
The electrical conductivity in his calculations was set 
at 1.0 x 10 5 mho/m resulting in a very large magnetic 
Reynolds number. In the present study, the calcu- 
lations are performed in the low magnetic Reynolds 
number regime using a realistic value of electrical con- 
ductivity. The flow is computed in two dimensions, 
but later will be extended to three dimensions. Be- 
cause of flow symmetry, only half of the channel is 
computed in the 2-D calculations. 

The flow in the nozzle section and the center sec- 
tion was computed using a combination of the OVER- 
FLOW code [26] and the present PNS code (without 
MHD effects). The initial conditions for the nozzle 
(flow at rest) were: 

po = 8.0 x 10 5 N/m 2 


T 0 = 7500 K 

The laminar flow was assumed to be in chemical equi- 
librium. The computed flowfield at the end of the 
center section was then used as the starting solu- 
tion for the flow calculation of the accelerator section. 
The MHD parameters used in the accelerator section 
were: 



= 50 mho/m 

By 

= 1.5 T 

E z 

= —Ku e B y 

i?e m 

= 0.05 


where the load factor (if) ranged in values from 0.0 to 
1.4, and the centerline velocity (u c ) at the beginning 
of the accelerator section had a value of 3162 m/s. 

The velocity profiles at the end of the accelerator 
section sure shown in Fig. 8 for different load factors. 
The velocity profile with no electric or magnetic fields 
is denoted by K = 0. The increase in the centerline 
velocity with distance (z) for various load factors is 
shown in Fig. 9. The centerline velocity increases 
by about 30% with a load factor of 1.4. It should 
be noted that the flow decelerates because of friction 
when no electric or magnetic fields are applied. 

Concluding Remarks 

In this study, a new parabolized Navier-Stokes al- 
gorithm has been developed to efficiently compute 
MHD flows in the low magnetic Reynolds number 
regime. The new algorithm has been used to com- 
pute both laminar and turbulent, supersonic, MHD 
flows over flat plates and in a rectangular accelera- 
tor section. Although only limited results have been 
obtained thus far, it can be seen that the present 
approach is quite promising. Computations of other 
test cases are currently underway in order to validate 
the current method. 
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Figure 6: Laminar/ turbulent skin friction coefficient 
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Figure 9: Centerline velocity for various load factors 
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